% -------------------------------------------------------------------------
% transformation from mee to eci
%
% [r, v] = mee2eci_(mu, mee)
%
% Input argumuents:
% -------------------------------------------------------------------------
% mu         [1x1 double]       gravitational constant    [ - ]        
% mee        [7x? double]       state variables           [ - ]
%
% Output argumuents:
% -------------------------------------------------------------------------
% r          [?x3 double]       position vector           [ - ]
% v          [?x3 double]       velocity vector           [ - ]
%
% External functions called:
% -------------------------------------------------------------------------
% None
%
% Copyright(C) 2015/06/22 by Chen Zhang, 
% School of Astronautics, Beihang University
% chenzhang.buaa@gmail.com
% -------------------------------------------------------------------------
function [r, v] = mee2eci_(mu, mee)
[~, col] = size(mee);
r = zeros(col , 3); v = zeros(col , 3);
p = mee(1 , :); ex = mee(2 , :); ey = mee(3 , :);
hx = mee(4 , :); hy = mee(5 , :); L = mee(6 , :);
cosl = cos(L);
sinl = sin(L);
W = 1 + ex .* cosl + ey .* sinl;
C = 1 + hx.^2 + hy.^2;
% compute eci position vector
r(: , 1) = p ./ W .* (cosl + (hx.^2 - hy.^2) .* ...
    cosl + 2 * hx .* hy .* sinl) ./ C;
r(: , 2) = p ./ W .* (sinl - (hx.^2 - hy.^2) .* ...
    sinl + 2 .* hx .* hy .* cosl) ./ C;
r(: , 3) = 2 * p ./ W .* (hx .* sinl - hy .* cosl) ./ C;
% compute eci velocity vector
v(: , 1) = - sqrt(mu ./ p) .* (sinl + (hx.^2 - hy.^2) .* sinl ...
    - 2 * hx .* hy .* cosl + ey - 2 * ex .* hx .* hy + ...
    (hx.^2 - hy.^2) .* ey) ./ C;
v(: , 2) = - sqrt(mu ./ p) .* (-cosl + (hx.^2 - hy.^2) .* cosl ...
    + 2 * hx .* hy .* sinl - ex + 2 * ey .* hx .* hy + ...
    (hx.^2 - hy.^2) .* ex) ./ C;
v(: , 3) = 2 * sqrt(mu ./ p) .* (hx .* cosl + hy .* sinl + ex .* hx ...
    + ey .* hy) ./ C;
end
% -------------------------------------------------------------------------
